# Replication material for:
# What drives public support for policies to enhance electric vehicle adoption?
# G. Brückmann & T. Benrnauer (2020) Environ. Res. Lett. 15 094002
# https://doi.org/10.1088/1748-9326/ab90a5

# Replication code by G. Brückmann
# Link to dataset: https://doi.org/10.7910/DVN/XF0GD9


if (!require("pacman")) install.packages("pacman")

pacman::p_load( 'cjoint',  'tidyverse', 'viridis',  'cregg', 'data.table')

set.seed(402)



sessionInfo()
# 
# R version 4.1.1 (2021-08-10)
# Platform: aarch64-apple-darwin20 (64-bit)
# Running under: macOS 13.3
# 
# Matrix products: default
# LAPACK: /Library/Frameworks/R.framework/Versions/4.1-arm64/Resources/lib/libRlapack.dylib
# 
# locale:
#   [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
# 
# attached base packages:
#   [1] grid      stats     graphics  grDevices utils     datasets  methods   base     
# 
# other attached packages:
#   [1] data.table_1.14.8 cregg_0.4.0       viridis_0.6.2     viridisLite_0.4.0 forcats_0.5.1     stringr_1.4.0     dplyr_1.1.1       purrr_0.3.4       readr_2.0.2       tidyr_1.1.4      
# [11] tibble_3.2.1      tidyverse_1.3.1   cjoint_2.1.0      survey_4.1-1      survival_3.2-11   Matrix_1.5-1      ggplot2_3.3.5     lmtest_0.9-38     zoo_1.8-9         sandwich_3.0-1   
# [21] pacman_0.5.1     
# 
# loaded via a namespace (and not attached):
#   [1] Rcpp_1.0.7       lubridate_1.7.10 lattice_0.20-44  assertthat_0.2.1 digest_0.6.28    utf8_1.2.2       mime_0.12        R6_2.5.1         cellranger_1.1.0 backports_1.2.1 
# [11] reprex_2.0.1     ggstance_0.3.5   httr_1.4.2       pillar_1.9.0     rlang_1.1.0      readxl_1.4.1     rstudioapi_0.13  splines_4.1.1    munsell_0.5.0    shiny_1.7.1     
# [21] broom_1.0.4      compiler_4.1.1   httpuv_1.6.3     modelr_0.1.8     pkgconfig_2.0.3  htmltools_0.5.5  mitools_2.4      tidyselect_1.2.0 gridExtra_2.3    fansi_0.5.0     
# [31] crayon_1.4.1     tzdb_0.1.2       dbplyr_2.1.1     withr_2.4.2      later_1.3.0      jsonlite_1.7.2   xtable_1.8-4     gtable_0.3.0     lifecycle_1.0.3  DBI_1.1.1       
# [41] magrittr_2.0.3   scales_1.2.1     cli_3.6.1        stringi_1.7.5    renv_0.17.3      fs_1.6.1         promises_1.2.0.1 xml2_1.3.3       ellipsis_0.3.2   generics_0.1.0  
# [51] vctrs_0.6.1      tools_4.1.1      glue_1.6.2       hms_1.1.1        fastmap_1.1.0    colorspace_2.0-2 rvest_1.0.1      haven_2.5.2 






################################ load data ################################ 
makeAllControls = FALSE

#Note: tt$evdt = 1 if the data comes from the EV sample dataset, tt$evdt==0 if the data comes from the non-EV sample. 

# loading itsself ===============
load("replicationdata_ERL2020GBTB.RData")
# in case the loaded data is not called tt anymore:
#tt <- replicationdata_ERL2020

# Make them all numeric:

for(i in names(tt)) {
  tt[[i]] = as.numeric(tt[[i]])
}

# Define the labels for the relevant variables
l_charging          = c("No new chargers", "At 1 of 1000 parkings", "At 10 of 1000 parkings", "At 100 of 1000 parkings")
l_subsidy       = c("No subsidy", "Subsidy CHF 1000", "Subsidy CHF 3000", "Subsidy CHF 5000")
l_registration  = c("New reg. allowed" , "New reg. forbidden")
l_info          = c("Keep current info level", "Increase info level", "Abandon infos")
l_financing   = c("No add. fin.","Savings","Increase income tax" ,"Malus for inefficient cars","Increase car vignette price")
l_hhinc         = c("Below CHF 4000","CHF 4000-8000","CHF 8000-12000","CHF 12000-16000","More than CHF 16000")
l_gender        = c("male","not male") 
l_droveev       = c("No","Yes") # tt$drovev:  : Non-EV Holder stated drove an EV before
l_age       = c("born < 1945","born 1945-1954", "born 1955-1964", "born 1965-1974", "born 1975-1984", "born >= 1985")
l_tesla = c("No","Yes") # tt$tesla : EV Holder with a tesla
l_icevbackup =c("No","Yes") # tt$icevbackup:  : EV Holder with a back-up non-bev car 





# Now put on the labels on the choice experimental variable
tt$charging     = factor(tt$charging,levels = c(0:3),labels = l_charging)
tt$subsidy      = factor(tt$subsidy,levels = c(0:3),labels = l_subsidy)
tt$registration = factor(tt$registration,levels=c(0,1),labels = l_registration)
tt$info         = factor(tt$info,levels=c(0,1,2),labels=l_info)
tt$financing    = factor(tt$financing,levels=c(0:4),labels = l_financing)

# Now the covariates
tt$gender          = factor(tt$gender,levels = c(1,2),labels = l_gender)
tt$hhinc        = factor(tt$hhinc,levels=c(1:5),labels=l_hhinc)
tt$everdrovebev = factor(tt$everdrovebev,levels = c(0,1),labels=l_droveev)
tt$Tesla = factor(tt$Tesla,levels = c(0,1),labels=l_tesla)
tt$icevbackup = factor(tt$icevbackup,levels = c(0,1),labels=l_icevbackup )
tt$age = factor(tt$age, levels=c(1:6),labels = l_age )



# Round the score from the environmental concern scale 
tt$envsc <- round(tt$envsc, digits = 0)

#subset if necessary: all those who did see information on financing visiblefin == 1
#tt_with_fin <- subset(tt, visiblefin == 1) 


# ====
# Missing value handling
# =====
# Mark complete cases (NOT individuals): which variables?
fullobs = c("choice","charging","subsidy","registration","info","rating")
tt$complete = complete.cases(tt[,fullobs,with=FALSE])


# Mark complete observations per ID and EV / ICEV dataset.
# If a person has a ten here, no missing values in all choice exp. observations
#tt  = tt[,c("full"):=.(sum(complete)),by=.(id,evdt)]

# Which one to keep?
tt = subset(tt,tt$complete==TRUE)

#####Creating a new unique identifier named UID

tt  =   tt[ , uid := .GRP, by = .(id,evdt)]



#pretty var.labels = c(info="Car information requirements", subsidy="EV purchase subsidy",
#               registration="Registration of highly fossil fuel consuming car",financing="Funding",
#              charging="Charging infrastructure provision at parkng space")


### the main data frame: 
df_conj <- tt
##### we will always reset to this.

## # ### useful functions: 
#GetLastAssigned und lastAssignedFnct 
GetLastAssigned <- function(match = "<- *data.frame",
                            remove = " *<-.*") {
  f <- tempfile()
  savehistory(f)
  history <- readLines(f)
  unlink(f)
  match <- grep(match, history, value = TRUE)
  get(sub(remove, "", match[length(match)]))
}

lastAssigned <- function(match = "<- *data.frame",
                         remove = " *<-.*") {
  f <- tempfile()
  savehistory(f)
  history <- readLines(f)
  unlink(f)
  match <- grep(match, history, value = TRUE)
  sub(remove, "", match[length(match)])
}

# what to plot tester
contentcheck<-function(x=(lastAssigned()))
{
  if(str_detect(((x)), "_c", negate = FALSE)) {
    is.combined <<- TRUE 
    is.ev <<- FALSE
  } else { is.combined <<- FALSE
  if(str_detect(((x)), "_a", negate = FALSE)) {
    is.pooled <<- TRUE 
    is.ev <<- FALSE
    is.combined <<- FALSE
  } else { is.pooled <<- FALSE
  if(str_detect(((x)), "_e", negate = FALSE)) {
    is.ev <<- TRUE 
    is.combined <<- FALSE
    is.pooled <<- FALSE
  } else {
    is.ev <<- FALSE
    is.combined <<- FALSE
    is.pooled <<- FALSE
  }
  }
  }
  if(str_detect(((x)), "_f", negate = FALSE)) {
    is.withfin <<- TRUE 
  } else {
    is.withfin  <<- FALSE 
  }  
  if(str_detect(((x)), "_t", negate = FALSE)) {
    is.withtes <<- TRUE 
  } else {
    is.withtes  <<- FALSE 
  }  
  if(str_detect(((x)), "_m", negate = FALSE)) {
    is.withmul <<- TRUE 
  } else {
    is.withmul  <<- FALSE 
  }  
}


#dataprepformula

prepare2<-function(x){
  require(data.table)
  x[,Group:=stringr::str_to_title(feature)]
  x[Group %chin% c('Info'), Group:=c('Information')]
  if(is.combined==TRUE){
    x[cars==1 ,Car:="No EV"]
    x[cars==2 ,Car:="EV"]
    oe <- c( 'No EV', 'EV', 'All')
    x[,"Car"]<-factor(x$Car, levels=oe)
    #   x <<- within(x, Car <- factor(Car, levels=oe))
  }else{
    if(is.ev==TRUE & is.pooled == FALSE){
      x[,Car:="EV"]
      oe <- c( 'No EV', 'EV', 'All')
      x[,"Car"] <- factor(x$Car, levels=oe)
      
    }else{
      if(is.ev==FALSE & is.pooled == FALSE){
        x[,Car:="No EV"]
        oe <- c( 'No EV', 'EV', 'All')
        x[,"Car"] <- factor(x$Car, levels=oe)
      }else {
        x[,Car:="All"]
        oe <- c( 'No EV', 'EV', 'All')
        x[,"Car"] <- factor(x$Car, levels=oe)
      }
    } 
  }
  if(is.withfin==TRUE){
    x[ ,Visible:="Yes"]
    visiblelable <- c("Yes", "No")
    x[,"Visible"] <- factor(x$Visible, levels=visiblelable)
    ogf <- c('Charging', 'Subsidy', 'Information', 'Registration', 'Financing' )
    olf <- rev(c('At 100 of 1000 parkings', 'At 10 of 1000 parkings', 'At 1 of 1000 parkings', 'No new chargers',
                 'Subsidy CHF 5000', 'Subsidy CHF 3000', 'Subsidy CHF 1000', 'No subsidy',
                 'Increase info level', 'Abandon infos', 'Keep current info level',
                 'New reg. forbidden', 'New reg. allowed',
                 'Increase car vignette price', 'Increase income tax', 'Malus for inefficient cars', 'Savings', 'No add. fin.'))
    x[,"Group"] <- factor(x$Group, levels=ogf)
    x[,"level"] <- factor(x$level, levels=olf)
  }else{
    x[ ,Visible:="No"]
    visiblelable <- c("Yes", "No")
    x[,"Visible"] <- factor(x$Visible, levels=visiblelable)
    ogf <- c('Charging', 'Subsidy', 'Information', 'Registration', 'Financing' )
    olf <- rev(c('At 100 of 1000 parkings', 'At 10 of 1000 parkings', 'At 1 of 1000 parkings', 'No new chargers',
                 'Subsidy CHF 5000', 'Subsidy CHF 3000', 'Subsidy CHF 1000', 'No subsidy',
                 'Increase info level', 'Abandon infos', 'Keep current info level',
                 'New reg. forbidden', 'New reg. allowed',
                 'Increase car vignette price', 'Increase income tax', 'Malus for inefficient cars', 'Savings', 'No add. fin.'))
    x[,"Group"] <- factor(x$Group, levels=ogf)
    x[,"level"] <- factor(x$level, levels=olf)
  } 
  return(x)
}

# for nice plotting, define s

set_s <- function(is.combined, is.withfin, is.ev, is.pooled){
  if(is.combined==FALSE & is.withfin==TRUE & is.ev==TRUE & is.pooled == FALSE){
    s <- 16
  } else if(is.combined==FALSE & is.withfin==TRUE & is.ev==FALSE & is.pooled == FALSE){
    s <- 17
  }
  else if(is.combined==FALSE & is.withfin==TRUE  & is.pooled == TRUE){
    s <- 15
  }
  else if(is.combined==FALSE & is.withfin==FALSE & is.pooled == TRUE){
    s <- 0
  }
  else if(is.combined==FALSE & is.withfin==FALSE & is.ev ==TRUE & is.pooled == FALSE){
    s <- 1
  }
  else if(is.combined==FALSE & is.withfin==FALSE & is.ev ==FALSE & is.pooled == FALSE){
    s <- 2
  }
  return(s)
}




# PRE ANALYSIS ------------------------------------------------------------------
#### Pre-Analysis
# check if cjoint finds the constraint in the data with funding visible: 
df_conj <- tt
setDT(df_conj)
df_conj <- df_conj[visiblefin==1]
# display frequencies
(f <- cj_freqs(df_conj, ~ charging +  subsidy +  registration + info + financing,
               id = ~uid))

# restrictions
## check display proportions
cj_props(df_conj, ~ financing, id = ~ uid)

## check which combinations were not allowed
subset(cj_props(df_conj, ~ subsidy + financing, id = ~ CaseID), Proportion == 0)
subset(cj_props(df_conj, ~ charging + financing, id = ~ CaseID), Proportion == 0)
# plotting
(p <- plot(f))
## change ggplot2 theme
p + ggplot2::theme_bw()

#Therefore, conclude that cjoint-package understood the restrictions. 








# ############ START OF THE ANALYSIS of all subgroups---------------------

# ##Let's first take a look at the main results (choice), later we will turn to rating.-------
# ######## SET ON CHOICE! ------------------------------------------------------------------
CHOICE= TRUE
RATING = FALSE
# Choice Financing Revealed ------------------------------------------------------------------
# Choice Financing Revealed A) ALL ------------------------------------------------------------------

df_conj <- tt
setDT(df_conj)
df_conj <- df_conj[visiblefin==1]
#this line gives us the relevant number of observations!
obsis_all_fin <- length(unique(df_conj$uid)) 
cmm_f_all <- data.frame(mm(df_conj, choice ~ charging +  subsidy +  registration + info + financing,
                           id = ~uid ))

lastAssigned()
contentcheck()

cmm_f_all<- setDT(cmm_f_all)
is.data.table(cmm_f_all)
cmm_f_all<- prepare2(cmm_f_all)


result <- set_s(is.combined, is.withfin, is.ev, is.pooled)
s <- print(result)

p_cmm_f_all      <- ggplot(cmm_f_all, aes(x = level, y = estimate, color=Car)) +
  geom_pointrange(aes(min = estimate - 1.95 * std.error, max = estimate + 1.95 * std.error), shape = s) +
  theme_bw() + 
  facet_wrap(Group ~ ., scales = "free_y", nrow = 5, strip.position = "left") +
  coord_flip() +
  ylab("Marginal Means") +
  xlab("Attributes") +
  ylim(0.2,0.7) +
  geom_hline(yintercept = 0.5, lty="dashed") +
  ggtitle(paste0("Binary Choice for Policy Bundles \nFinancing Revealed (All car holders, n = ", obsis_all_fin,")")) +
  theme(plot.title = element_text(size = 11, face = "bold"))

p_cmm_f_all
g <- ggplotGrob(p_cmm_f_all)
g$layout$l[g$layout$name == "title"] <- 1
p_cmm_f_all  <- grid::grid.draw(g)

p_cmm_f_all


# Choice Financing Revealed B) ICEV ------------------------------------------------------------------

df_conj <- tt
setDT(df_conj)
df_conj <- df_conj[visiblefin==1&EV==0]

#this line gives us the relevant number of observations!
obsis <- length(unique(df_conj$uid)) 

cmm_f_icev <- data.frame(mm(df_conj, choice ~ charging +  subsidy +  registration + info + financing,
                            id = ~uid ))


lastAssigned()
contentcheck()

setDT(cmm_f_icev)
cmm_f_icev <- prepare2(cmm_f_icev)

result <- set_s(is.combined, is.withfin, is.ev, is.pooled)
s <- print(result)

p_cmm_f_icev      <- ggplot(cmm_f_icev, aes(x = level, y = estimate, color=Car)) +
  geom_pointrange(aes(min = estimate - 1.95 * std.error, max = estimate + 1.95 * std.error), shape = s) +
  theme_bw() + 
  facet_wrap(Group ~ ., scales = "free_y", nrow = 5, strip.position = "left") +
  coord_flip() +
  ylab("Marginal Means") +
  xlab("Attributes") +
  ylim(0.2,0.7) +
  geom_hline(yintercept = 0.5, lty="dashed") +
  ggtitle(paste0("Binary Choice for Policy Bundles \nFinancing Revealed (No-EV holders, n =", obsis,")")) +
  #  geom_text(aes(label = Car), colour = "black", size = 2.5, hjust=1.05, vjust=1.2)
  if(is.pooled==FALSE){scale_colour_manual(values=c( "#330042","#4EC150"),aesthetics = "colour")}

p_cmm_f_icev



# Choice Financing Revealed C) EV ------------------------------------------------------------------

df_conj <- tt
setDT(df_conj)
df_conj <- df_conj[visiblefin==1&EV==1]

#this line gives us the relevant number of observations!
obsis <- length(unique(df_conj$uid)) 

cmm_f_ev <- data.frame(mm(df_conj, choice ~ charging +  subsidy +  registration + info + financing,
                          id = ~uid ))

lastAssigned()
contentcheck()

setDT(cmm_f_ev)
cmm_f_ev <- prepare2(cmm_f_ev)

result <- set_s(is.combined, is.withfin, is.ev, is.pooled)
s <- print(result)


p_cmm_f_ev  <- ggplot(cmm_f_ev, aes(x = level, y = estimate, color = Car)) +
  geom_pointrange(aes(min = estimate - 1.95 * std.error, max = estimate + 1.95 * std.error), shape = s) +
  theme_bw() + 
  facet_wrap(Group ~ ., scales = "free_y", nrow = 5, strip.position = "left") +
  coord_flip() +
  ylab("Marginal Means") +
  xlab("Attributes") +
  ylim(0.2,0.7) +
  scale_colour_viridis_d(begin = 0.75 , end = 0) +
  geom_hline(yintercept = 0.5, lty="dashed") +
  ggtitle(paste0("Binary Choice for Policy Bundles \n Financing Revealed (EV holders, n = ",obsis,")")) 

p_cmm_f_ev


# Choice Financing Revealed d) ICEV & EV ------------------------------------------------------------------
## FIGURE 5 in paper
cmm_f_c <- data.frame(bind_rows(cmm_f_icev, cmm_f_ev, .id = "cars"))

lastAssigned()
contentcheck()

setDT(cmm_f_c)
cmm_f_c<- prepare2(cmm_f_c)


p_cmm_f_c      <- ggplot(cmm_f_c, aes(x = level, y = estimate, color=Car)) +
  geom_pointrange(aes(min = estimate - 1.95 * std.error, max = estimate + 1.95 * std.error, shape = Car)) +
  theme_bw() + 
  facet_wrap(Group ~ ., scales = "free_y", nrow = 5, strip.position = "left") +
  coord_flip() +
  ylab("Marginal Means") +
  xlab("Attributes") +
  ylim(0.2,0.7) +
  geom_hline(yintercept = 0.5, lty="dashed") +
  ggtitle(paste0("Binary Choice for Policy Bundles \n Financing Revealed (All car holders, n =", obsis_all_fin,")")) 
if(is.combined==TRUE & is.withfin == TRUE){ p_cmm_f_c <- p_cmm_f_c + scale_shape_manual(values=c(17,16))} 
if(is.pooled==FALSE){p_cmm_f_c <- p_cmm_f_c + scale_colour_manual(values=c( "#330042","#4EC150"), aesthetics = "colour")} 
if(is.combined==TRUE & is.withfin==FALSE){p_cmm_f_c <- p_cmm_f_c +  scale_shape_manual(values=c(2,1))}
#  geom_text(aes(label = Car), colour = "black", size = 2.5, hjust=1.05, vjust=1.2)
geom_text(aes(label = Car), colour = "black", size = 2.5, hjust=1.05, vjust=1.2)
p_cmm_f_c
# # End of FIGURE 5 in paper #



# 2)  fin hidden -----------------------------------------------------
# Choice Financing Hidden A) ALL ------------------------------------------------------------------

# resetting:
df_conj <- tt
setDT(df_conj)
df_conj <- df_conj[visiblefin==0]

#this line gives us the relevant number of observations!
obsis_all_hidden <- length(unique(df_conj$uid)) 

cmm_h_all <- data.frame(mm(df_conj, choice ~ charging +  subsidy +  registration + info ,
                           id = ~uid ))


lastAssigned()
contentcheck()

cmm_h_all<- setDT(cmm_h_all)
is.data.table(cmm_h_all)
cmm_h_all<- prepare2(cmm_h_all)


result <- set_s(is.combined, is.withfin, is.ev, is.pooled)
s <- print(result)

p_cmm_h_all      <- ggplot(cmm_h_all, aes(x = level, y = estimate, color=Car)) +
  geom_pointrange(aes(min = estimate - 1.95 * std.error, max = estimate + 1.95 * std.error), shape = s) +
  theme_bw() + 
  facet_wrap(Group ~ ., scales = "free_y", nrow = 5, strip.position = "left") +
  coord_flip() +
  ylab("Marginal Means") +
  xlab("Attributes") +
  ylim(0.2,0.7) +
  geom_hline(yintercept = 0.5, lty="dashed") +
  ggtitle(paste0("Binary Choice for Policy Bundles \nFinancing Hidden (All car holders, n =", obsis_all_hidden ,")")) +
  #  geom_text(aes(label = Car), colour = "black", size = 2.5, hjust=1.05, vjust=1.2)
  if(is.pooled==FALSE){scale_colour_manual(col="#4EC150", col="#330042")} 
p_cmm_h_all


# Choice Financing Hidden B) ICEV ------------------------------------------------------------------
# # FIGURE 2 in paper
df_conj <- tt
setDT(df_conj)
df_conj <- df_conj[visiblefin==0&EV==0]

#this line gives us the relevant number of observations!
obsis <- length(unique(df_conj$uid)) 

cmm_h_icev <- data.frame(mm(df_conj, choice ~ charging +  subsidy +  registration + info ,
                            id = ~uid ))

lastAssigned()
contentcheck()

setDT(cmm_h_icev)
cmm_h_icev <- prepare2(cmm_h_icev)

result <- set_s(is.combined, is.withfin, is.ev, is.pooled)
s <- print(result)

p_cmm_h_icev    <- ggplot(cmm_h_icev, aes(x = level, y = estimate, color=Car)) +
  geom_pointrange(aes(min = estimate - 1.95 * std.error, max = estimate + 1.95 * std.error), shape = s) +
  theme_bw() + 
  facet_wrap(Group ~ ., scales = "free_y", nrow = 5, strip.position = "left") +
  coord_flip() +
  ylab("Marginal Means") +
  xlab("Attributes") +
  ylim(0.2,0.7) +
  geom_hline(yintercept = 0.5, lty="dashed") +
  ggtitle(paste0("Binary Choice for Policy Bundles \nFinancing Revealed (No-EV holders, n =",obsis,")" )) +
  #  geom_text(aes(label = Car), colour = "black", size = 2.5, hjust=1.05, vjust=1.2)
  if(is.pooled==FALSE){scale_colour_manual(values=c( "#330042","#4EC150"),aesthetics = "colour")}

p_cmm_h_icev


####nicely
result <- set_s(is.combined, is.withfin, is.ev, is.pooled)
s <- print(result)

p_cmm_h_icev_n    <- ggplot(cmm_h_icev, aes(x = level, y = estimate, color=Car)) +
  geom_pointrange(aes(min = estimate - 1.95 * std.error, max = estimate + 1.95 * std.error), shape = s) +
  theme_bw() + 
  facet_wrap(Group ~ ., scales = "free_y", nrow = 5, strip.position = "left") +
  coord_flip() +
  ylab("Marginal Means") +
  xlab("Attributes") +
  ylim(0.2,0.7) +
  geom_hline(yintercept = 0.5, lty="dashed") +
  ggtitle(paste0("Binary Choice for Policy Bundles \nFinancing Revealed (No-EV holders, n =",obsis,")" )) +
  guides(colour="none", size="none", shape = "none") +
  theme(legend.position="bottom")+
  if(is.pooled==FALSE){scale_colour_manual(values=c( "#330042","#4EC150"),aesthetics = "colour")}

p_cmm_h_icev_n
# # End of FIGURE 2 in paper #




# Choice Financing Hidden C) EV ------------------------------------------------------------------

df_conj <- tt
setDT(df_conj)
df_conj <- df_conj[visiblefin==0&EV==1]

#this line gives us the relevant number of observations!
obsis <- length(unique(df_conj$uid)) 

cmm_h_ev <- data.frame(mm(df_conj, choice ~ charging +  subsidy +  registration + info ,
                          id = ~uid ))


lastAssigned()
contentcheck()

setDT(cmm_h_ev)
cmm_h_ev <- prepare2(cmm_h_ev)

result <- set_s(is.combined, is.withfin, is.ev, is.pooled)
s <- print(result)


p_cmm_h_ev    <- ggplot(cmm_h_ev, aes(x = level, y = estimate, color=Car)) +
  geom_pointrange(aes(min = estimate - 1.95 * std.error, max = estimate + 1.95 * std.error), shape = s) +
  theme_bw() + 
  facet_wrap(Group ~ ., scales = "free_y", nrow = 5, strip.position = "left") +
  coord_flip() +
  ylab("Marginal Means") +
  xlab("Attributes") +
  ylim(0.2,0.7) +
  geom_hline(yintercept = 0.5, lty="dashed") +
  ggtitle(paste0("Binary Choice for Policy Bundles \nFinancing Hidden (EV holders, n =", obsis, ")" )) +
  scale_colour_viridis_d(begin = 0.75 , end = 0) 

p_cmm_h_ev




# Choice Financing Revealed d) ICEV & EV ------------------------------------------------------------------
# FIGURE 4 in paper #
cmm_h_c <- data.frame(bind_rows(cmm_h_icev, cmm_h_ev, .id = "cars"))

lastAssigned()
contentcheck()

setDT(cmm_h_c)
cmm_h_c<- prepare2(cmm_h_c)

p_cmm_h_c      <- ggplot(cmm_h_c, aes(x = level, y = estimate, color=Car)) +
  geom_pointrange(aes(min = estimate - 1.95 * std.error, max = estimate + 1.95 * std.error, shape = Car)) +
  theme_bw() + 
  facet_wrap(Group ~ ., scales = "free_y", nrow = 5, strip.position = "left") +
  coord_flip() +
  ylab("Marginal Means") +
  xlab("Attributes") +
  ylim(0.2,0.7) +
  geom_hline(yintercept = 0.5, lty="dashed") +
  ggtitle(paste0("Binary Choice for Policy Bundles \nFinancing Revealed (All car holders, n =", obsis_all_hidden,")")) 
if(is.combined == TRUE & is.withfin == TRUE){ p_cmm_h_c <- p_cmm_h_c + scale_shape_manual(values=c(17,16))} 
if(is.pooled == FALSE){p_cmm_h_c <- p_cmm_h_c + scale_colour_manual(values=c( "#330042","#4EC150"),aesthetics = "colour")} 
if(is.combined == TRUE & is.withfin == FALSE){p_cmm_h_c <- p_cmm_h_c + scale_shape_manual(values=c(2,1))}

p_cmm_h_c







##### For Figure 3 in Paper
cmm_c <- data.frame(bind_rows( cmm_f_icev, cmm_h_icev, .id = "funding"))


setDT(cmm_c)

lastAssigned()
contentcheck()

setDT(cmm_c)

preparecarcombi<-function(x){
  require(data.table)
  x[,Group:=stringr::str_to_title(feature)]
  x[Group %chin% c('Info'), Group:=c('Information')]
  if(is.combined==TRUE){
    x[funding==1 ,Visible:="Yes"]
    x[funding==2 ,Visible:="No"]
    oe <- c( 'No EV', 'EV', 'All')
    x[,"Car"]<-factor(x$Car, levels=oe)
    visiblelable <- c("Yes", "No")
    x[,"Visible"] <- factor(x$Visible, levels=visiblelable)
    ogf <- c('Charging', 'Subsidy', 'Information', 'Registration', 'Financing' )
    olf <- rev(c('At 100 of 1000 parkings', 'At 10 of 1000 parkings', 'At 1 of 1000 parkings', 'No new chargers',
                 'Subsidy CHF 5000', 'Subsidy CHF 3000', 'Subsidy CHF 1000', 'No subsidy',
                 'Increase info level', 'Abandon infos', 'Keep current info level',
                 'New reg. forbidden', 'New reg. allowed',
                 'Increase car vignette price', 'Increase income tax', 'Malus for inefficient cars', 'Savings', 'No add. fin.'))
    x[,"Group"] <- factor(x$Group, levels=ogf)
    x <- x[order(ogf),]
    x[,"level"] <- factor(x$level, levels=olf)
  }else{
    print("ups!")}
  return(x) }

preparecarcombi(cmm_c)




p_cmm_f_c      <- ggplot(cmm_c, aes(x = level, y = estimate, color=Car, shape = Visible)) +
  geom_linerange(aes(min = estimate - 1.95 * std.error, max = estimate + 1.95 * std.error), size=.6) +
  geom_point(aes(size=Visible)) +
  theme_bw() + 
  facet_wrap(Group ~ ., scales = "free_y", nrow = 5, strip.position = "left") +
  coord_flip() +
  ylab("Marginal Means") +
  xlab("Attributes") +
  ylim(0.2,0.7)+
  scale_colour_manual(values=c( "#330042","#4EC150"),aesthetics = "colour")+
  geom_hline(yintercept = 0.5, lty="dashed") +
  ggtitle("Binary Choice for Policy Bundles (No-EV holders, n = 3768)") +
  scale_shape_manual(values=c(17,2))+
  scale_size_manual(values=c(3,3)) +
  labs(shape="Financing Visible") +
  guides(shape=guide_legend(override.aes=list(size=3, linetype=0)),
         size="none", colour = "none")+
  theme(legend.position="bottom")

p_cmm_f_c
# # End of FIGURE 3 in paper #

















# #### FOR THE APPENDIX: ------------------------------------------------------------------


# Set to Rating ------------------------------------------------------------------

RATING= TRUE
CHOICE = FALSE


# Rating Financing Revealed ------------------------------------------------------------------
df_conj <- tt
setDT(df_conj)
df_conj <- df_conj[visiblefin==1]

#this line gives us the relevant number of observations!
obis_all_fin <- length(unique(df_conj$uid)) 

# Rating Financing Revealed A) ALL ------------------------------------------------------------------

rmm_f_all <- data.frame(mm(df_conj, rating ~ charging +  subsidy +  registration + info + financing,
                           id = ~uid ))




lastAssigned()
contentcheck()

rmm_f_all<- setDT(rmm_f_all)
is.data.table(rmm_f_all)
rmm_f_all<- prepare2(rmm_f_all)

result <- set_s(is.combined, is.withfin, is.ev, is.pooled)
s <- print(result)

p_rmm_f_all      <- ggplot(rmm_f_all, aes(x = level, y = estimate, color=Car)) +
  geom_pointrange(aes(min = estimate - 1.95 * std.error, max = estimate + 1.95 * std.error), shape = s) +
  theme_bw() + 
  facet_wrap(Group ~ ., scales = "free_y", nrow = 5, strip.position = "left") +
  coord_flip() +
  ylab("Marginal Means") +
  xlab("Attributes") +
  ggtitle("Rating of Policy Bundles \n Financing Revealed (All car holders, n = 2442)") +
  if(is.pooled==FALSE){scale_colour_manual(col="#4EC150", col="#330042")} 
p_rmm_f_all


# Rating Financing Revealed B) ICEV ------------------------------------------------------------------
df_conj <- tt
setDT(df_conj)
df_conj <- df_conj[visiblefin==1&EV==0]

#this line gives us the relevant number of observations!
obsis <- length(unique(df_conj$uid)) 

rmm_f_icev <- data.frame(mm(df_conj, rating ~ charging +  subsidy +  registration + info + financing,
                            id = ~uid ))

#plot(rmm_f_icev)

lastAssigned()
contentcheck()

setDT(rmm_f_icev)
rmm_f_icev <- prepare2(rmm_f_icev)


result <- set_s(is.combined, is.withfin, is.ev, is.pooled)
s <- print(result)

p_rmm_f_icev      <- ggplot(rmm_f_icev, aes(x = level, y = estimate, color=Car)) +
  geom_pointrange(aes(min = estimate - 1.95 * std.error, max = estimate + 1.95 * std.error), shape = s) +
  theme_bw() + 
  facet_wrap(Group ~ ., scales = "free_y", nrow = 5, strip.position = "left") +
  coord_flip() +
  ylab("Marginal Means") +
  xlab("Attributes") +
  ggtitle("Rating of Policy Bundles \n Financing Revealed (No-EV holders, n = 1871)") +
  if(is.pooled==FALSE){scale_colour_manual(values=c( "#330042","#4EC150"),aesthetics = "colour")}

p_rmm_f_icev


# Rating Financing Revealed C) EV ------------------------------------------------------------------

df_conj <- tt
setDT(df_conj)
df_conj <- df_conj[visiblefin==1&EV==1]

#this line gives us the relevant number of observations!
obsis <- length(unique(df_conj$uid)) 


rmm_f_ev <- data.frame(mm(df_conj, rating ~ charging +  subsidy +  registration + info + financing,
                          id = ~uid ))

#plot(rmm_f_ev)

lastAssigned()
contentcheck()

setDT(rmm_f_ev)
rmm_f_ev <- prepare2(rmm_f_ev)

result <- set_s(is.combined, is.withfin, is.ev, is.pooled)
s <- print(result)

p_rmm_f_ev  <- ggplot(rmm_f_ev, aes(x = level, y = estimate, color = Car)) +
  geom_pointrange(aes(min = estimate - 1.95 * std.error, max = estimate + 1.95 * std.error), shape = s) +
  theme_bw() + 
  facet_wrap(Group ~ ., scales = "free_y", nrow = 5, strip.position = "left") +
  coord_flip() +
  ylab("Marginal Means") +
  xlab("Attributes") +
  #  ylim(-0.3,0.3)+
  scale_colour_viridis_d(begin = 0.75 , end = 0) +
  ggtitle("Rating of Policy Bundles \n Financing Revealed (EV holders, n = 571)") 
p_rmm_f_ev



# Rating Financing Revealed d) ICEV & EV ------------------------------------------------------------------
### BEGINNING OF SI FIG 6
rmm_f_c <- data.frame(bind_rows(rmm_f_icev, rmm_f_ev, .id = "cars"))

lastAssigned()
contentcheck()

setDT(rmm_f_c)
rmm_f_c<- prepare2(rmm_f_c)


p_rmm_f_c      <- ggplot(rmm_f_c, aes(x = level, y = estimate, color=Car)) +
  geom_pointrange(aes(min = estimate - 1.95 * std.error, max = estimate + 1.95 * std.error, shape = Car)) +
  theme_bw() + 
  facet_wrap(Group ~ ., scales = "free_y", nrow = 5, strip.position = "left") +
  coord_flip() +
  ylab("Marginal Means") +
  xlab("Attributes") +
  ggtitle(paste0("Rating of Policy Bundles \n Financing Revealed (All car holders, n =", obsis_all_fin,")")) 
if(is.combined==TRUE & is.withfin == TRUE){ p_rmm_f_c <- p_rmm_f_c + scale_shape_manual(values=c(17,16))} 
if(is.pooled==FALSE){p_rmm_f_c <- p_rmm_f_c + scale_colour_manual(values=c( "#330042","#4EC150"),aesthetics = "colour")} 
if(is.combined==TRUE & is.withfin==FALSE){p_rmm_f_c <- p_rmm_f_c +  scale_shape_manual(values=c(2,1))}

p_rmm_f_c

### END OF SI FIG 6

# 2)  fin hidden -----------------------------------------------------
# Rating Financing Hidden A) ALL ------------------------------------------------------------------
# resetting:

df_conj <- tt
setDT(df_conj)
df_conj <- df_conj[visiblefin==0]

#this line gives us the relevant number of observations!
obsis_all_nofin <- length(unique(df_conj$uid)) 

rmm_h_all <- data.frame(mm(df_conj, rating ~ charging +  subsidy +  registration + info ,
                           id = ~uid ))


lastAssigned()
contentcheck()

rmm_h_all<- setDT(rmm_h_all)
is.data.table(rmm_h_all)
rmm_h_all<- prepare2(rmm_h_all)

result <- set_s(is.combined, is.withfin, is.ev, is.pooled)
s <- print(result)

p_rmm_h_all      <- ggplot(rmm_h_all, aes(x = level, y = estimate, color=Car)) +
  geom_pointrange(aes(min = estimate - 1.95 * std.error, max = estimate + 1.95 * std.error), shape = s) +
  theme_bw() + 
  facet_wrap(Group ~ ., scales = "free_y", nrow = 5, strip.position = "left") +
  coord_flip() +
  ylab("Marginal Means") +
  xlab("Attributes") +
  ggtitle(paste0("Rating of Policy Bundles \n Financing Hidden (All car holders, n =", obsis_all_nofin,")")) +
  if(is.pooled==FALSE){scale_colour_manual(col="#4EC150", col="#330042")} 
p_rmm_h_all



# Rating Financing Hidden B) ICEV ------------------------------------------------------------------
#### Beginning FIG SI 3
df_conj <- tt
setDT(df_conj)
df_conj <- df_conj[visiblefin==0&EV==0]

#this line gives us the relevant number of observations!
obsis <- length(unique(df_conj$uid)) 

rmm_h_icev <- data.frame(mm(df_conj, rating ~ charging +  subsidy +  registration + info ,
                            id = ~uid ))

lastAssigned()
contentcheck()

setDT(rmm_h_icev)
rmm_h_icev <- prepare2(rmm_h_icev)

result <- set_s(is.combined, is.withfin, is.ev, is.pooled)
s <- print(result)

p_rmm_h_icev    <- ggplot(rmm_h_icev, aes(x = level, y = estimate, color=Car)) +
  geom_pointrange(aes(min = estimate - 1.95 * std.error, max = estimate + 1.95 * std.error), shape = s) +
  theme_bw() + 
  facet_wrap(Group ~ ., scales = "free_y", nrow = 5, strip.position = "left") +
  coord_flip() +
  ylab("Marginal Means") +
  xlab("Attributes") +
  ylim(2,5)+
  ggtitle(paste0("Rating of Policy Bundles \n Financing Revealed (No-EV holders, n =", obsis, ")")) +
  if(is.pooled==FALSE){scale_colour_manual(values=c( "#330042","#4EC150"),aesthetics = "colour")}

p_rmm_h_icev
#### END FIG SI 3

# Rating Financing Hidden C) EV ------------------------------------------------------------------

df_conj <- tt
setDT(df_conj)
df_conj <- df_conj[visiblefin==0&EV==1]

#this line gives us the relevant number of observations!
obsis <- length(unique(df_conj$uid)) 

rmm_h_ev <- data.frame(mm(df_conj, rating ~ charging +  subsidy +  registration + info ,
                          id = ~uid ))


lastAssigned()
contentcheck()

setDT(rmm_h_ev)
rmm_h_ev <- prepare2(rmm_h_ev)

result <- set_s(is.combined, is.withfin, is.ev, is.pooled)
s <- print(result)

p_rmm_h_ev    <- ggplot(rmm_h_ev, aes(x = level, y = estimate, color=Car)) +
  geom_pointrange(aes(min = estimate - 1.95 * std.error, max = estimate + 1.95 * std.error), shape = s) +
  theme_bw() + 
  facet_wrap(Group ~ ., scales = "free_y", nrow = 5, strip.position = "left") +
  coord_flip() +
  ylab("Marginal Means") +
  xlab("Attributes") +
  ggtitle("Rating of Policy Bundles \n Financing Revealed (EV holders, n = 571)") +
  if(is.pooled==FALSE){scale_colour_manual(values=c( "#4EC150","#330042"),aesthetics = "colour")}

p_rmm_h_ev



# Rating Financing Revealed d) ICEV & EV ------------------------------------------------------------------
### FIG SI 5:
rmm_h_c <- data.frame(bind_rows(rmm_h_icev, rmm_h_ev, .id = "cars"))

lastAssigned()
contentcheck()

setDT(rmm_h_c)
rmm_h_c<- prepare2(rmm_h_c)


p_rmm_h_c      <- ggplot(rmm_h_c, aes(x = level, y = estimate, color=Car)) +
  geom_pointrange(aes(min = estimate - 1.95 * std.error, max = estimate + 1.95 * std.error, shape = Car)) +
  theme_bw() + 
  facet_wrap(Group ~ ., scales = "free_y", nrow = 5, strip.position = "left") +
  coord_flip() +
  ylab("Marginal Means") +
  xlab("Attributes") +
  ggtitle(paste0("Rating of Policy Bundles \n Financing Hidden (All car holders, n =", obsis_all_hidden,")")) 
if(is.combined == TRUE & is.withfin == TRUE){ p_rmm_h_c <- p_rmm_h_c + scale_shape_manual(values=c(17,16))} 
if(is.pooled == FALSE){p_rmm_h_c <- p_rmm_h_c + scale_colour_manual(values=c( "#330042","#4EC150"),aesthetics = "colour")} 
if(is.combined == TRUE & is.withfin == FALSE){p_rmm_h_c <- p_rmm_h_c + scale_shape_manual(values=c(2,1))}

p_rmm_h_c
### END FIG SI 5


##### For SI Figure 4 in Paper
rmm_c <- data.frame(bind_rows( rmm_f_icev, rmm_h_icev, .id = "funding"))


setDT(rmm_c)

lastAssigned()
contentcheck()

setDT(rmm_c)


preparecarcombi(rmm_c)




p_rmm_f_c      <- ggplot(rmm_c, aes(x = level, y = estimate, color=Car, shape = Visible)) +
  geom_linerange(aes(min = estimate - 1.95 * std.error, max = estimate + 1.95 * std.error), size=.6) +
  geom_point(aes(size=Visible)) +
  theme_bw() + 
  facet_wrap(Group ~ ., scales = "free_y", nrow = 5, strip.position = "left") +
  coord_flip() +
  ylab("Marginal Means") +
  xlab("Attributes") +
 # ylim(0.2,0.7)+
  scale_colour_manual(values=c( "#330042","#4EC150"),aesthetics = "colour")+
  ggtitle("Rating of Policy Bundles (No-EV holders, n = 3768)") +
  scale_shape_manual(values=c(17,2))+
  scale_size_manual(values=c(3,3)) +
  labs(shape="Financing Visible") +
  guides(shape=guide_legend(override.aes=list(size=3, linetype=0)),
         size="none", colour = "none")+
  theme(legend.position="bottom")

p_rmm_f_c
# # End of FIGURE SI 4 in paper #





#MORE APPENDIX


# 1) with fin revealed -----------------------------------------------------
df_conj <- tt
setDT(df_conj)
df_conj <- df_conj[visiblefin==1]

hyp_df_wf = expand.grid(charging = unique(df_conj$charging), subsidy = unique(df_conj$subsidy),
                        registration = unique(df_conj$registration), info = unique(df_conj$info), 
                        financing = unique(df_conj$financing),  EV = unique(df_conj$EV)
)

#----------
out_m_wf = lm(rating  ~ (charging + subsidy + registration + info + financing) * EV, data = df_conj)

hyp_df_wf$pred_out = predict(out_m_wf, newdata = hyp_df_wf)

hyp_df_wf$visiblefin =1


#### FOR SI TABLE 1 FIND IN THIS DF THE HIGHEST VALUE FOR THOSE WHO SAW FUNDING
highest <- hyp_df_wf[hyp_df_wf$EV==0,]


fig_rate_average_wf <- ggplot(hyp_df_wf, aes(x=pred_out, colour=factor(EV, labels = c("No-EV","EV") ) ) )+ 
  geom_density(size=1.2) + 
  theme_bw() +
  theme(panel.grid.major.x = element_line(size=.5, color="gray80"),
        panel.grid.minor.y = element_blank(),
        panel.grid.minor.x = element_blank(),
        panel.grid.major.y = element_line(size=.5, color="gray80" )) +
  xlab("Distribution of average ratings for \n policy-packages financing revealed") +
  ylab("Density") +
  scale_color_manual(values = c("#330042","#4EC150")) +
  scale_x_continuous(breaks = c(1,2,3,4,5,6,7), labels = c("1","2", "3", "4", "5", "6", "7")) +
  geom_vline(xintercept = 4, lty="dashed") +
  geom_vline(xintercept = 5, lty="dashed")+
  labs(colour="Car")

fig_rate_average_wf



# 2)  fin hidden -----------------------------------------------------

# resetting:
df_conj <- tt
setDT(df_conj)
df_conj <- df_conj[visiblefin==0]

hyp_df_fh = expand.grid(charging = unique(df_conj$charging), subsidy = unique(df_conj$subsidy),
                        registration = unique(df_conj$registration), info = unique(df_conj$info), 
                        EV = unique(df_conj$EV)
)

#----------
out_m_fh = lm(rating  ~ (charging + subsidy + registration + info ) * EV, data = df_conj)

hyp_df_fh$pred_out = predict(out_m_fh, newdata = hyp_df_fh)
hyp_df_fh$visiblefin = 0


#### FOR SI TABLE 1 FIND IN THIS DF THE HIGHEST VALUE FOR THOSE WHO DID NOT SAW FUNDING
highest_hidden <- hyp_df_fh[hyp_df_fh$EV==0,]


fig_rate_average_fh <- ggplot(hyp_df_fh, aes(x=pred_out, colour=factor(EV, labels = c("No-EV","EV") ) ) )+
  geom_density(size=1.2) + 
  theme_bw() +
  theme(panel.grid.major.x = element_line(size=.5, color="gray80"),
        panel.grid.minor.y = element_blank(),
        panel.grid.minor.x = element_blank(),
        panel.grid.major.y = element_line(size=.5, color="gray80" )) +
  xlab("Distribution of average ratings for \n policy-packages financing revealed") +
  ylab("Density") +
  scale_color_manual(values = c("#330042","#4EC150")) +
  scale_x_continuous(breaks = c(1,2,3,4,5,6,7), labels = c("1","2", "3", "4", "5", "6", "7")) +
  geom_vline(xintercept = 4, lty="dashed") +
  geom_vline(xintercept = 5, lty="dashed") +
  labs(colour="Car")
fig_rate_average_fh


# 3)  Join the 2 figures  -----------------------------------------------------
#### END OF APPENDIX FIG 2
setDT(hyp_df_fh)
setDT(hyp_df_wf)
hyp_df_all <- merge(hyp_df_wf, hyp_df_fh, all = TRUE)
setDT(hyp_df_all)



fig_rate_average_all_scaled <- ggplot(hyp_df_all, aes(x=pred_out,  ..scaled.., colour=factor(interaction(EV,visiblefin), labels = c("No-EV, No", "EV, No", "No-EV, Yes","EV, Yes" )) )) + 
  geom_density(size=1.2) + 
  theme_bw() +
  theme(panel.grid.major.x = element_line(size=.5, color="gray80"),
        panel.grid.minor.y = element_blank(),
        panel.grid.minor.x = element_blank(),
        panel.grid.major.y = element_line(size=.5, color="gray80" )) +
  xlab("Average rating (1-7 Likert-scale)") +
  ylab("Density") +
  scale_colour_manual(values = c("No-EV, No" = "#330042", "EV, No" = "#4EC150", "No-EV, Yes" = "#26547B","EV, Yes" = "#FCE51E"  )) +
  scale_x_continuous(breaks = c(1,2,3,4,5,6,7), labels = c("1","2", "3", "4", "5", "6", "7")) +
  geom_vline(xintercept = 4, lty="dashed") +
  geom_vline(xintercept = 5, lty="dashed") +
  labs(colour="Car Type and Financing Visible") +
  ggtitle("Distribution of predicted average ratings for policy-packages")

fig_rate_average_all_scaled
# ggsave("Average ratings by car and financing scaled.pdf",fig_rate_average_all_scaled,width = 160,height=180,units = "mm",scale = 1.2)
#### END OF APPENDIX FIG 2

# 3)  And now for how many packages pass certain thresholds  -----------------------------------------------------
dat.m <- melt(hyp_df_all,id.vars=c("pred_out", "EV", "visiblefin"), measure.vars=c("charging", "subsidy", "registration", "info", "financing"))
hyp_df <- hyp_df_all
nrow(hyp_df[hyp_df$pred_out >= 4 & hyp_df$EV == 0 & hyp_df$visiblefin == 0, ])
nrow(hyp_df[hyp_df$pred_out >= 4 & hyp_df$EV == 0 & hyp_df$visiblefin == 1,])
nrow(hyp_df[hyp_df$pred_out >= 4 & hyp_df$EV == 1 & hyp_df$visiblefin == 0,])
nrow(hyp_df[hyp_df$pred_out >= 4 & hyp_df$EV == 1 & hyp_df$visiblefin == 1,])

nrow(hyp_df[hyp_df$pred_out >= 5 & hyp_df$EV == 0 & hyp_df$visiblefin == 0, ])
nrow(hyp_df[hyp_df$pred_out >= 5 & hyp_df$EV == 0 & hyp_df$visiblefin == 1 ,])
nrow(hyp_df[hyp_df$pred_out >= 5 & hyp_df$EV == 1 & hyp_df$visiblefin == 0 ,])
nrow(hyp_df[hyp_df$pred_out >= 5 & hyp_df$EV == 1 & hyp_df$visiblefin == 1,])

# very unliked packages
nrow(hyp_df[hyp_df$pred_out <= 2.5 & hyp_df$EV == 0 & hyp_df$visiblefin == 0, ])
nrow(hyp_df[hyp_df$pred_out <= 2.5 & hyp_df$EV == 0 & hyp_df$visiblefin == 1 ,])
nrow(hyp_df[hyp_df$pred_out <= 2.5 & hyp_df$EV == 1 & hyp_df$visiblefin == 0 ,])
nrow(hyp_df[hyp_df$pred_out <=2.5 & hyp_df$EV == 1 & hyp_df$visiblefin == 1,])

nrow(hyp_df[hyp_df$pred_out <= 3 & hyp_df$EV == 0 & hyp_df$visiblefin == 0, ])
nrow(hyp_df[hyp_df$pred_out <= 3 & hyp_df$EV == 0 & hyp_df$visiblefin == 1 ,])
nrow(hyp_df[hyp_df$pred_out <= 3 & hyp_df$EV == 1 & hyp_df$visiblefin == 0 ,])
nrow(hyp_df[hyp_df$pred_out <= 3 & hyp_df$EV == 1 & hyp_df$visiblefin == 1,])

nrow(hyp_df[hyp_df$pred_out <= 3.1 & hyp_df$EV == 0 & hyp_df$visiblefin == 0, ])
nrow(hyp_df[hyp_df$pred_out <= 3.1 & hyp_df$EV == 0 & hyp_df$visiblefin == 1 ,])
nrow(hyp_df[hyp_df$pred_out <= 3.1 & hyp_df$EV == 1 & hyp_df$visiblefin == 0 ,])
nrow(hyp_df[hyp_df$pred_out <= 3.1 & hyp_df$EV == 1 & hyp_df$visiblefin == 1,])

### Until here part of appendix
